from numpy import array,linspace,zeros,log
from matplotlib.pyplot import plot,subplot,title,show
from numpy.random import random
from numpy.linalg import norm

def Fx(x):
    return 1;

def Fy(y):
    return y;

#x0=0;
#y0=1;
t0=0;
tFinal=1;
T=(10**4)*5;
t=linspace(t0,tFinal,T)
dt=t[1]-t[0]
#print(dt)
#x=zeros(T)
#y=zeros(T)
#x[0]=x0;
#y[0]=y0;

#for i in range(1,T):
#    x[i]=x[i-1]+Fx(x[i-1])*dt;
#    y[i]=y[i-1]+Fy(y[i-1])*dt;
    
#erroX=abs(x-t);
#erroY=abs(y-exp(t))
#print(erroY)

#print(dt)
#subplot(2,3,1)
#title('x(t)')
#plot(t,x,'r')
#subplot(2,3,2)
#title('y(t)')
#plot(t,y,'b')
#subplot(2,3,4)
#title('erro x(t)')
#plot(t,erroX,'r')
#subplot(2,3,5)
#title('erro y(t)')
#plot(t,erroY,'b')

J=array([0,0,0,2]).reshape(2,2);
norma=[]
delta=random(2)*(10**-10)
#delta0=random.random(2)*(10**-6)
#delta=zeros(2)
#delta[0]=delta0[0];
#delta[1]=delta0[1];
norma.append(norm(delta))
norma[0]=1;
delta=delta/norm(delta)
y2=zeros(T)
Y=zeros(T)

fracao=0
normalizacao=1
for i in range(1,T):
    delta=delta+dt*J*delta
    #print(delta)    
    norma.append(norm(delta))
    y2[i]=((i-1)/(i))*y2[i-1]+2*log(norma[i]/norma[i-1])
    if i>fracao*T:
        Y[i]=((i-1)*Y[i-1]+y2[i])/i;
    
    if i%normalizacao==0:
        delta=delta/norm(delta)
        norma[i]=1;
        
    #delta=delta/norm(delta)
    norma[i]=1;
    
subplot(2,2,1)    
title('y')
plot(t,y2)
subplot(2,2,2)
title('Y')
plot(t,Y)
subplot(2,2,3)
title('erro y')
plot(t,abs(y2-2*t))
subplot(2,2,4)
title('erro Y')
plot(t,abs(Y-t))
#plot(t,abs(Y-t))
show()